Sumário Executivo
A Ibema, juntamente com a Microsoft, solicitaram a BlueShift que avaliasse a possibilidade de melhorar as estratégias de go-to-market da Ibema por meio da aplicação de técnicas de inteligência artificial e de machine learning aos dados de volumes de vendas e de preços da Ibema.
A fim de atender a essa demanda, a BlueShift criou modelos de forecasting de volumes e de preços líquidos para a Ibema por meio de técnicas tanto consagradas quanto de ponta, a fim de apoiar o processo decisório da empresa no que diz respeito à otimização de esforços e de alocação de recursos para regiões e sites de maior interesse, ou seja, aqueles com melhor prognóstico de faturamento.
A aplicação das técnicas mencionadas anteriormente para os dados do estado de São Paulo, bobinas do site TVO, permite concluir que os modelos preditivos desenvolvidos apresentaram muito bom desempenho e que, portanto, podem ser usados para prever acuradamente tanto o preço líquido quanto o volume de vendas em um horizonte de até 14 dias. Essas previsões, se estendidas para outras regiões, sites e tipos de embalagem, certamente poderão apoiar o processo decisório da Ibema no que diz respeito ao planejamento e controle da produção, bem como em relação às estratégias de marketing.
Introdução
O objetivo desta PoC é criar modelos de forecasting de volumes e de preços líquidos para a Ibema por meio de técnicas tanto consagradas quanto de ponta, a fim de apoiar o processo decisório da empresa no que diz respeito à otimização de esforços e de recursos para regiões e sites de maior interesse.
Será usada como métrica de interesse de uma região/site a previsão de faturamento líquido - preço líquido x volume - nos próximos 14 dias.
Ingestão e Preparação dos Dados
Os dados recebidos da Ibema em formato xlsx foram ingeridos e, depois de alguns esforços de análise exploratória, receberam o seguinte tratamento:
Variáveis tidas por categóricas:
cod_site,cod_produto,cod_uf,cod_tipo_embalagem,cod_gramatura,cod_configuracaoVariáveis tomadas como estritamente positivas:
num_comprimento,num_largura,qtd_volume,vlr_faturamento_bruto,vlr_preco_liquido
Variáveis não utilizadas: vlr_despesa_financeira, vlr_faturamento_bruto, vlr_comissao, vlr_despesa_venda, vlr_frete_unitario, nom_cidade, vlr_faturamento_sem_imposto, vlr_impostos, vlr_taxa_nota, vlr_faturamento_sem_ipi, vlr_frete
A seguir, mostram-se algumas estatísticas descritivas do dataset resultante dessa etapa, considerando:
- UF: ‘SP’
- SITE: ‘TVO’
- EMBALAGEM: ‘BOBINA’
UF <- 'SP'
SITE <- 'TVO'
EMBALAGEM <- 'BOBINA'
ibema <- read_delim(
"00_data/Dados_BlueShift_Total.csv", ";",
escape_double = FALSE,
col_types = cols(
VLR_PRECO_LIQUIDO = col_number(),
VLR_FATURAMENTO_BRUTO = col_number(),
VLR_FATURAMENTO_SEM_IPI = col_number(),
VLR_FATURAMENTO_SEM_IMPOSTO = col_number(),
VLR_TAXA_NOTA = col_number(),
QTD_VOLUME = col_number()
),
trim_ws = TRUE
) %>%
janitor::clean_names() %>%
mutate(
cod_gramatura = as.factor(cod_gramatura),
cod_site = as.factor(cod_site),
cod_uf = as.factor(cod_uf),
cod_tipo_embalagem = as.factor(cod_tipo_embalagem),
cod_configuracao = as.factor(cod_configuracao),
cod_produto = as.factor(cod_produto),
dat_movimento = as.Date(dat_movimento),
cod_produto = fct_lump_prop(cod_produto, prop = 0.01),
cod_configuracao = fct_lump_prop(cod_configuracao, prop = 0.01),
cod_uf = fct_lump_prop(cod_uf, prop = 0.01),
cod_gramatura = fct_lump_prop(cod_gramatura, prop = 0.01)
) %>%
filter(
num_comprimento > 0, num_largura > 0, qtd_volume > 0,
vlr_faturamento_bruto > 0, vlr_preco_liquido > 0
) %>%
select(
-vlr_despesa_financeira, -vlr_faturamento_bruto,
-vlr_comissao, -vlr_despesa_venda, -vlr_frete_unitario,
-nom_cidade, -vlr_faturamento_sem_imposto, -vlr_impostos,
-vlr_taxa_nota, -vlr_faturamento_sem_ipi, -vlr_frete
) %>%
filter(cod_uf == UF, cod_site == SITE, cod_tipo_embalagem == EMBALAGEM) %>%
select(-cod_uf, -cod_site, -cod_tipo_embalagem) %>%
arrange(dat_movimento) %>%
pad_by_time(.date_var = dat_movimento, .by = "day", .pad_value = 0) %>%
summarise_by_time(.date_var = dat_movimento, .by = "day",
preco_liq_medio = mean(vlr_preco_liquido),
volume_total = sum(qtd_volume)) %>%
mutate(
preco_liq_medio = zoo::na.locf(preco_liq_medio),
volume_total = replace_na(volume_total, 0)
)
skim_without_charts(ibema)| Name | ibema |
| Number of rows | 1179 |
| Number of columns | 3 |
| _______________________ | |
| Column type frequency: | |
| Date | 1 |
| numeric | 2 |
| ________________________ | |
| Group variables | None |
Variable type: Date
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| dat_movimento | 0 | 1 | 2018-01-03 | 2021-03-26 | 2019-08-15 | 1179 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 |
|---|---|---|---|---|---|---|---|---|---|
| preco_liq_medio | 0 | 1 | 2364.38 | 1951.63 | 0 | 0 | 3469.03 | 4027.81 | 5168.57 |
| volume_total | 0 | 1 | 23.37 | 28.98 | 0 | 0 | 13.48 | 38.13 | 195.79 |
Feature Engineering
A etapa anterior resultou em uma série temporal cujos últimos 365 dias são mostrados a seguir, a título de ilustração:
ibema %>%
tail(365) %>%
plot_time_series(.date_var = dat_movimento, .value = volume_total)Da análise visual dessa série temporal, algumas características se destacam, a saber:
Parece haver uma sazonalidade semanal
Há picos não sazonais cuja explicação pode demandar o cruzamento desses dados com outros dados, por exemplo: eventos de marketing
Para confrontar as hipóteses descritas anteriormente, mostram-se a seguir os gráficos de autocorrelação e de autocorrelação parcial:
ibema %>%
plot_acf_diagnostics(.date_var = dat_movimento,
.value = log1p(volume_total))ibema %>%
plot_seasonal_diagnostics(.date_var = dat_movimento,
.value = log1p(volume_total))Da análise dos gráficos anteriores, nota-se a sazonalidade semanal (7 dias) e também um certo padrão sazonal ao longo dos dias da semana. Essas informações serão importantes no processo de construção dos modelos matemáticos.
Forecasting
calibrate_and_plot <- function(..., type = "testing") {
if (type == "testing") {
new_data <- testing(splits)
} else {
new_data <- training(splits) %>% drop_na()
}
calibration_tbl <- modeltime_table(...) %>%
modeltime_calibrate(new_data)
print(modeltime_accuracy(calibration_tbl))
calibration_tbl %>%
modeltime_forecast(
new_data = new_data,
actual_data = data_prepared_tbl
) %>%
plot_modeltime_forecast(.conf_interval_show = FALSE)
}
refit_and_plot <- function(...) {
calibration_tbl <- modeltime_table(...) %>%
modeltime_calibrate(testing(splits))
refit_tbl <- calibration_tbl %>%
modeltime_refit(data_prepared_tbl)
refit_tbl %>%
modeltime_forecast(new_data = forecast_tbl,
actual_data = data_prepared_tbl) %>%
# Invert Transformation
mutate(across(.value:.conf_hi, .fns = ~ standardize_inv_vec(
x = .,
mean = std_mean,
sd = std_sd
))) %>%
mutate(across(.value:.conf_hi, .fns = ~ log_interval_inv_vec(
x = .,
limit_lower = limit_lower,
limit_upper = limit_upper,
offset = offset
))) %>%
plot_modeltime_forecast()
}Transformação dos dados
Em face das análises anteriores, podem ser usados os lags de 7, 14 e 21 dias como tentativas de variáveis adicionais, além do dia da semana e do mês, visando melhorar o desempenho dos modelos matemáticos.
Adicionalmente, a fim de melhorar o desempenho dos modelos matemáticos de machine learning, tomou-se o logaritmo do volume diário, que posteriormente foi normalizado (média zero e desvio padrão igual a um).
# Transformação dos dados
ibema_transformed_tbl <- ibema %>%
# Preprocessamento do Alvo
mutate(vol_tot_trans = log_interval_vec(volume_total, limit_lower = 0,
offset = 1)) %>%
mutate(vol_tot_trans = standardize_vec(vol_tot_trans)) %>%
select(-volume_total, -preco_liq_medio)
# Parâmetros chave
limit_lower <- 0
limit_upper <- 216.3745
offset <- 1
std_mean <- -3.17308718389508
std_sd <- 1.92403817890689
ibema_transformed_tbl %>% glimpse()## Rows: 1,179
## Columns: 2
## $ dat_movimento <date> 2018-01-03, 2018-01-04, 2018-01-05, 2018-01-06, 2018-01~
## $ vol_tot_trans <dbl> 0.77029164, -1.14306006, -1.14306006, -1.14306006, -1.14~
Visualização do efeito de feature engineering
A fim de testar as hipóteses que motivaram as transformações de dados descritas anteriormente, a seguir mostra-se uma regressão linear do volume (transformado) em função das datas e de:
lags: 7, 14 e 21 dias
dia da semana: de segunda-feira a domingo
mês: de janeiro a dezembro
plot_time_series_regression(
.data = ibema_transformed_tbl %>%
tk_augment_lags(.value = vol_tot_trans, .lags = c(7, 14, 21)) %>%
drop_na(),
.date_var = dat_movimento,
.formula = vol_tot_trans ~ as.numeric(dat_movimento)
+ wday(dat_movimento, label = TRUE)
+ month(dat_movimento, label = TRUE)
+ .
- dat_movimento,
.show_summary = TRUE
)##
## Call:
## stats::lm(formula = .formula, data = .data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.76844 -0.31520 0.03258 0.48072 2.37820
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.866e-01 1.279e+00 0.224 0.822748
## as.numeric(dat_movimento) -1.593e-05 7.052e-05 -0.226 0.821363
## wday(dat_movimento, label = TRUE).L 2.354e-01 6.210e-02 3.792 0.000158
## wday(dat_movimento, label = TRUE).Q -1.602e+00 1.004e-01 -15.954 < 2e-16
## wday(dat_movimento, label = TRUE).C 8.591e-02 6.105e-02 1.407 0.159663
## wday(dat_movimento, label = TRUE)^4 -5.819e-01 6.759e-02 -8.609 < 2e-16
## wday(dat_movimento, label = TRUE)^5 -2.290e-03 6.088e-02 -0.038 0.970007
## wday(dat_movimento, label = TRUE)^6 -6.125e-02 6.089e-02 -1.006 0.314676
## month(dat_movimento, label = TRUE).L 7.952e-02 8.026e-02 0.991 0.321998
## month(dat_movimento, label = TRUE).Q 8.877e-02 8.179e-02 1.085 0.277968
## month(dat_movimento, label = TRUE).C -8.202e-02 8.079e-02 -1.015 0.310188
## month(dat_movimento, label = TRUE)^4 -2.399e-01 8.010e-02 -2.995 0.002802
## month(dat_movimento, label = TRUE)^5 8.740e-02 8.012e-02 1.091 0.275564
## month(dat_movimento, label = TRUE)^6 -4.255e-02 8.018e-02 -0.531 0.595767
## month(dat_movimento, label = TRUE)^7 -1.950e-02 7.912e-02 -0.247 0.805337
## month(dat_movimento, label = TRUE)^8 3.009e-04 7.927e-02 0.004 0.996972
## month(dat_movimento, label = TRUE)^9 2.785e-02 8.097e-02 0.344 0.730918
## month(dat_movimento, label = TRUE)^10 -1.385e-01 8.148e-02 -1.700 0.089418
## month(dat_movimento, label = TRUE)^11 1.390e-01 8.184e-02 1.698 0.089788
## vol_tot_trans_lag7 2.316e-02 2.986e-02 0.776 0.438083
## vol_tot_trans_lag14 -3.405e-02 2.989e-02 -1.139 0.254863
## vol_tot_trans_lag21 -3.311e-02 2.977e-02 -1.112 0.266311
##
## (Intercept)
## as.numeric(dat_movimento)
## wday(dat_movimento, label = TRUE).L ***
## wday(dat_movimento, label = TRUE).Q ***
## wday(dat_movimento, label = TRUE).C
## wday(dat_movimento, label = TRUE)^4 ***
## wday(dat_movimento, label = TRUE)^5
## wday(dat_movimento, label = TRUE)^6
## month(dat_movimento, label = TRUE).L
## month(dat_movimento, label = TRUE).Q
## month(dat_movimento, label = TRUE).C
## month(dat_movimento, label = TRUE)^4 **
## month(dat_movimento, label = TRUE)^5
## month(dat_movimento, label = TRUE)^6
## month(dat_movimento, label = TRUE)^7
## month(dat_movimento, label = TRUE)^8
## month(dat_movimento, label = TRUE)^9
## month(dat_movimento, label = TRUE)^10 .
## month(dat_movimento, label = TRUE)^11 .
## vol_tot_trans_lag7
## vol_tot_trans_lag14
## vol_tot_trans_lag21
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7831 on 1136 degrees of freedom
## Multiple R-squared: 0.3994, Adjusted R-squared: 0.3883
## F-statistic: 35.97 on 21 and 1136 DF, p-value: < 2.2e-16
Do sumário e do gráfico anteriores, depreende-se que os trabalhos de feature engineering resultaram em um excelente ponto de partida para os algoritmos de machine learning.
Métricas de comparação de desempenho
Daqui em diante, serão testadas nove abordagens para prever os volumes dos próximos 14 dias e os modelos resultantes serão comparados, principalmente, pelas seguintes métricas:
MAE: Mean Absolute Error
RMSE: Root Mean Squared Error
RSQ: R2
Resultados
O gráfico a seguir mostra os dados de treinamento e os de teste (14 dias) usados para treinar e verificar o desempenho dos modelos de forecasting selecionados.
# 1.0 PASSO 1 - CRIAR FULL DATA SET
# - Estender pela janela futura
# - Feature engineering
# - Regressores externos
horizon <- 14 # 14 dias
lag_period <- 14
rolling_periods <- c(7, 14, 21)
dot_value <- str_glue("vol_tot_trans_lag{lag_period}")
data_prepared_full_tbl <- ibema_transformed_tbl %>%
# Janela futura
bind_rows(
future_frame(.data = ., .date_var = dat_movimento, .length_out = horizon)
) %>%
# Lags autocorrelacionados
tk_augment_lags(vol_tot_trans, .lags = lag_period) %>%
# Rolling features to stabilize trend from M5 Competition
tk_augment_slidify(
.value = rlang::sym(dot_value),
.f = mean,
.period = rolling_periods,
.align = "center",
.partial = TRUE
) %>%
# Format Columns
rename_with(.cols = contains("lag"), .fn = ~ str_c("lag_", .x))
# 2.0 PASSO 2 - SEPARAR os DADOS DE MODELAGEM E DE FORECAST
data_prepared_tbl <- data_prepared_full_tbl %>%
filter(!is.na(vol_tot_trans))
forecast_tbl <- data_prepared_full_tbl %>%
filter(is.na(vol_tot_trans))
# 3.0 TRAIN/TEST (DATASET DE MODELAGEM)
splits <- time_series_split(data_prepared_tbl, assess = horizon,
cumulative = TRUE)
splits %>%
tk_time_series_cv_plan() %>%
plot_time_series_cv_plan(dat_movimento, vol_tot_trans)# 4.0 RECIPES
# - Time Series Signature - Adds bulk time-based features
# - More Transformation to Time Series Signature vars
# - Interactions, if any
# - Fourier Features
recipe_spec_base <- recipe(vol_tot_trans ~ ., data = training(splits)) %>%
# Time Series Signature
step_timeseries_signature(dat_movimento) %>%
step_rm(matches("(iso)|(xts)|(hour)|(minute)|(second)|(am.pm)")) %>%
step_rm(ends_with("(day)|(day7)")) %>%
# Standardization
step_normalize(matches("(index.num)|(year)")) %>%
# Dummy Encoding (One Hot Encoding)
step_dummy(all_nominal(), one_hot = TRUE) %>%
# Fourier
step_fourier(dat_movimento, period = rolling_periods, K = 2)
# recipe_spec_base %>% prep() %>% juice() %>% glimpse()
# 5.0 LINEAR REGRESSION MODEL
# * Model spec
model_spec_lm <- linear_reg() %>%
set_engine("lm")
# * Model Workflow
recipe_spec_lm <- recipe_spec_base %>%
step_rm(dat_movimento) %>%
step_naomit(contains("lag"))
workflow_fit_lm <- workflow() %>%
add_model(model_spec_lm) %>%
add_recipe(recipe_spec_lm) %>%
fit(training(splits))
# 6.0 GLM REGRESSION MODEL
# * GLM Model Spec
model_spec_glm <- linear_reg(penalty = 0.1, mixture = 0.5) %>%
set_engine("glmnet")
# * GLM Model Workflow
recipe_spec_glm <- recipe_spec_base %>%
step_rm(dat_movimento) %>%
step_naomit(contains("lag"))
workflow_fit_glm <- workflow() %>%
add_model(model_spec_glm) %>%
add_recipe(recipe_spec_glm) %>%
fit(training(splits))
# 7.0 RANDOM FOREST MODEL
# * RF Model Spec
model_spec_rf <- rand_forest() %>%
set_engine("randomForest")
# * RF Model Workflow
workflow_fit_rf <- workflow() %>%
add_model(model_spec_rf) %>%
add_recipe(recipe_spec_glm) %>%
fit(training(splits))
# 8.0 ARIMA MODEL
model_fit_auto_arima <- arima_reg() %>%
set_engine("auto_arima") %>%
fit(
vol_tot_trans ~ dat_movimento
+ fourier_vec(dat_movimento, period = 7)
+ fourier_vec(dat_movimento, period = 14)
+ fourier_vec(dat_movimento, period = 21),
data = training(splits)
)
# 9.0 PROPHET MODEL, NO XREGS
model_fit_prophet <- prophet_reg(
# changepoint_num = 25,
# changepoint_range = 0.8,
# seasonality_yearly = FALSE,
# seasonality_weekly = TRUE,
# seasonality_daily = TRUE
) %>%
set_engine("prophet") %>%
fit(vol_tot_trans ~ dat_movimento, data = training(splits))
# 10.0 PROPHET MODEL, XREGS
model_fit_prophet_xregs <- prophet_reg(
# changepoint_num = 25,
# changepoint_range = 0.8,
# seasonality_yearly = FALSE,
# seasonality_weekly = TRUE,
# seasonality_daily = TRUE
) %>%
set_engine("prophet") %>%
fit(
vol_tot_trans ~ dat_movimento
+ fourier_vec(dat_movimento, period = 7)
+ fourier_vec(dat_movimento, period = 14)
+ fourier_vec(dat_movimento, period = 21),
data = training(splits)
)
# 11.0 CUBIST
model_spec_cubist <- cubist_rules() %>%
set_engine("Cubist")
set.seed(123)
workflow_fit_cubist <- workflow() %>%
add_model(model_spec_cubist) %>%
add_recipe(recipe_spec_lm) %>%
fit(training(splits))
# 12.0 XGBOOST
workflow_fit_xgboost <- workflow() %>%
add_model(
boost_tree() %>%
set_engine("xgboost")
) %>%
add_recipe(recipe_spec_lm) %>%
fit(training(splits))
# 13.0 NNETAR
model_spec_nnetar <- nnetar_reg(
# non_seasonal_ar = 1,
# seasonal_ar = 3,
# hidden_units = 10,
# penalty = 10,
# num_networks = 20,
# epochs = 50
) %>%
set_engine("nnetar")
set.seed(123)
workflow_fit_nnetar <- workflow() %>%
add_model(model_spec_nnetar) %>%
add_recipe(recipe_spec_base) %>%
fit(training(splits) %>% drop_na())A tabela a seguir mostra os resultados do desempenho dos modelos matemáticos treinados nos 14 dias de teste, ou seja, nos 14 dias que não participaram do treinamento dos algoritmos.
Os algoritmos comparados a seguir são os seguintes:
Regressão linear
GLMNET: modelos lineares generalizados via máxima verossimilhança penalizada
Random Forest
ARIMA: média móvel autoregressiva integrada
PROPHET: algoritmo desenvolvido e utilizado pelo Facebook
CUBIST: algoritmo baseado em árvores, como o Random Forest, mas com equações lineares nos nós terminais
XGBOOST: Extreme Gradient Boosting
NNETAR: uma combinação de redes neurais e de modelos autorregressivos
# 13.0 MODELTIME WORKFLOW
# * Modeltime Table
submodels_tbl <- modeltime_table(
workflow_fit_lm, workflow_fit_glm, workflow_fit_rf,
model_fit_auto_arima, model_fit_prophet,
model_fit_prophet_xregs, workflow_fit_cubist,
workflow_fit_xgboost, workflow_fit_nnetar
)
# * Calibrate Testing Data
submodels_calibrated_tbl <- submodels_tbl %>%
modeltime_calibrate(testing(splits))
# * Measure Test Accuracy
submodels_calibrated_tbl %>%
modeltime_accuracy() %>%
table_modeltime_accuracy(
.interactive = TRUE,
bordered = TRUE,
resizable = TRUE
)Modelo Ensemble
Da tabela anterior, excetuando-se o modelo ARIMA, notam-se os excelentes desempenhos dos demais modelos. Uma estratégia para construir um modelo ainda melhor, a partir dos treinados até o momento, consiste em criar um modelo ensemble que calcule suas previsões a partir da média dos três “melhores” modelos, a saber: RANDOMFOREST, XGBOOST e PROPHET.
Os resultados dessa estratégia são mostrados a seguir:
# 15.0 ENSEMBLE
# * Make Ensemble
ensemble_fit_mean <- submodels_tbl %>%
filter(.model_desc %in% c("RANDOMFOREST", "XGBOOST", "PROPHET")) %>%
ensemble_average(type = "mean")
# * Modeltime Table
ensemble_tbl <- modeltime_table(
ensemble_fit_mean
)
# * Ensemble Test Accuracy
ensemble_tbl %>%
combine_modeltime_tables(
submodels_tbl %>%
filter(.model_desc %in% c("RANDOMFOREST", "XGBOOST", "PROPHET"))
) %>%
modeltime_accuracy(testing(splits)) %>%
table_modeltime_accuracy(
.interactive = TRUE,
bordered = TRUE,
resizable = TRUE
)Finalmente, os gráficos a seguir mostram, respectivamente, as previsões de volumes nos próximos 14 dias calculadas pelo modelo ensemble, comparadas com os dados de teste, bem como nos próximos 14 dias do horizonte de previsão:
ensemble_tbl %>%
modeltime_calibrate(testing(splits)) %>%
modeltime_forecast(
new_data = testing(splits),
actual_data = data_prepared_tbl,
keep_data = TRUE
) %>%
# Invert Transformation
mutate(across(.value:.conf_hi, .fns = ~ standardize_inv_vec(
x = .,
mean = std_mean,
sd = std_sd
))) %>%
mutate(across(.value:.conf_hi, .fns = ~ log_interval_inv_vec(
x = .,
limit_lower = limit_lower,
limit_upper = limit_upper,
offset = offset
))) %>%
plot_modeltime_forecast()ensemble_refit_tbl <- ensemble_tbl %>%
modeltime_calibrate(testing(splits)) %>%
modeltime_refit(data_prepared_tbl)
# * Visualize Ensemble Forecast
ensemble_refit_tbl %>%
modeltime_forecast(
new_data = forecast_tbl,
actual_data = data_prepared_tbl,
keep_data = TRUE
) %>%
# Invert Transformation
mutate(across(.value:.conf_hi, .fns = ~ standardize_inv_vec(
x = .,
mean = std_mean,
sd = std_sd
))) %>%
mutate(across(.value:.conf_hi, .fns = ~ log_interval_inv_vec(
x = .,
limit_lower = limit_lower,
limit_upper = limit_upper,
offset = offset
))) -> fore_vol
plot_modeltime_forecast(fore_vol)Forecasting dos Preços Líquidos
Como o objetivo desta PoC é prever o faturamento líquido nos próximos 14 dias - em um estado, site e para uma determinada embalagem - , nos faltam as previsões correspondentes dos preços líquidos.
A bem da concisão, não repetiremos todo o processo anterior para os preços líquidos, cujos dados passaram por processo semelhante ao apresentado neste relatório para os volumes diários.
Modelou-se a evolução dos preços líquidos com as mesmas técnicas e algoritmos utilizados para os volumes, cujos resultados constam na tabela a seguir:
precos <- read_rds("00_artifacts/price_artifacts.rds")
# * Calibrate Testing Data ----
submodels_calibrated_tbl <- precos$submodels_tbl %>%
modeltime_calibrate(precos$test_tbl)
# * Measure Test Accuracy ----
submodels_calibrated_tbl %>%
modeltime_accuracy() %>%
table_modeltime_accuracy(
.interactive = TRUE,
bordered = TRUE,
resizable = TRUE
)Da análise da tabela anterior, depreende-se que o melhor modelo é o NNETAR, sendo que os gráficos a seguir mostram, respectivamente, as previsões de preços líquidos nos próximos 14 dias calculadas pelo modelo selecionado, comparadas com os dados de teste, bem como nos próximos 14 dias do horizonte de previsão:
best_model_fit <- precos$best_model_fit
calibration_tbl <- modeltime_table(best_model_fit) %>%
modeltime_calibrate(new_data = precos$test_tbl)
calibration_tbl %>%
modeltime_forecast(
new_data = precos$test_tbl,
actual_data = precos$data_prepared_tbl
) %>%
# Invert Transformation
mutate(across(.value:.conf_hi, .fns = ~ standardize_inv_vec(
x = .,
mean = precos$inv_params$std_mean,
sd = precos$inv_params$std_sd
))) %>%
mutate(across(.value:.conf_hi, .fns = exp
)) %>%
plot_modeltime_forecast()refit_tbl <- calibration_tbl %>%
modeltime_refit(precos$data_prepared_tbl)
# Visualize Forecast
refit_tbl %>%
modeltime_forecast(
new_data = precos$fore_tbl,
actual_data = precos$data_prepared_tbl,
keep_data = TRUE
) %>%
# Invert Transformation
mutate(across(.value:.conf_hi, .fns = ~ standardize_inv_vec(
x = .,
mean = precos$inv_params$std_mean,
sd = precos$inv_params$std_sd
))) %>%
mutate(across(.value:.conf_hi, .fns = exp
)) -> fore_price
plot_modeltime_forecast(fore_price)Forecasting do Faturamento Líquido
De posse das previsões de volumes e de preços, calculam-se as previsões de faturamento por meio da multiplicação de ambas, conforme mostrado no gráfico a seguir:
fore_vol %>%
tail(180) %>%
# dplyr::filter(.key == "prediction") %>%
select(dat_movimento, .key, .value) %>%
rename(vol = .value) %>%
bind_cols(
fore_price %>%
tail(180) %>%
# dplyr::filter(.key == "prediction") %>%
select(.value) %>%
rename(price = .value)
) %>%
mutate(net_income = vol * price) %>%
pivot_longer(cols = c(vol, price, net_income)) %>%
plot_time_series(.date_var = dat_movimento, .value = value,
.color_var = .key, .facet_vars = name, .smooth = FALSE)Conclusões e Recomendações
A aplicação das técnicas mencionadas anteriormente para os dados do estado de São Paulo, bobinas do site TVO, permite concluir que os modelos preditivos desenvolvidos apresentaram muito bom desempenho e que, portanto, podem ser usados para prever acuradamente tanto o preço líquido quanto o volume de vendas em um horizonte de até 14 dias. Essas previsões, se estendidas para outras regiões, sites e tipos de embalagem, certamente poderão apoiar o processo decisório da Ibema no que diz respeito ao planejamento e controle da produção, bem como em relação às estratégias de marketing.
Em face das conclusões anteriores, recomenda-se o seguinte:
Avaliar modelos para outras unidades da federação, tipos de embalagem, cidades, etc.
Listar os forecastings por ordem decrescente de faturamento líquido previsto, conforme a granularidade que se deseje avaliar, e.g. estado, cidade, tipo de embalagem, configuração de produto, etc.
Automatizar o processo por meio de um sistema HPFS - High Performance Forecasting System